R Markdown
if (!require("knitr")) install.packages("knitr")
if (!require("SummarizedExperiment")) install.packages("SummarizedExperiment")
if (!require("edgeR")) install.packages("edgeR")
if (!require("tidyr")) install.packages("tidyr")
if (!require("scales")) install.packages("scales")
if (!require("plotly")) install.packages("plotly")
# IMPORT VARIABLES
data_samples <- params$data_samples
data_counts <- params$data_counts
data_annotation <- params$data_annotation
setGeneName <- params$setGeneName
cpm_value <- params$cpm_value
excluded_samples <- params$excluded_samples
design_value <- params$design_value
matrix_value <- params$matrix_value
alpha <- params$alpha
# Get desgin
design_value <- reOrderDesign(matrix_value, design_value, data_samples)
# FILTER DATA IF NESSECARY
#Filter if nessecary
data_samples <- data_samples[!rownames(data_samples) %in% excluded_samples,]
data_counts <- data_counts[,!colnames(data_counts) %in% excluded_samples]
# SHOW VALUES OF VARIABLES
cpm_value
## [1] 1
excluded_samples
## [1] ""
design_value
## [1] "~0 + phenotype"
matrix_value
## [1] "pDC - pDC_unstimulated"
alpha
## [1] 0.05
# SHOW RAW DATA
se <- readCountsFromTable(data_counts, data_samples)
alignmentSummaryPlot(se)
complexityPlot(se)
# READ ALL FILES INTO SE
se <- readCountsFromTable(data_counts[!grepl('^__', rownames(data_counts)),], data_samples)
se <- addSamplesFromTableToSE(se, data_samples)
if (!is.null(data_annotation)) {
se <- addAnnotationsFromTableToSE(se, data_annotation)
}
# GET AND CREATE DGE
dge <- DGEList(counts = assay(se), samples = colData(se), genes = rowData(se))
row.names(dge$genes) <- row.names(dge$counts)
dge <- dge[ rowSums( abs( dge$counts ) ) > 1, ]
tempDge <- dge
tempDge$counts <- cpm(dge, log = TRUE)
countDistributionLinePlot(tempDge)
# GET SELECTED FEATURES
edger <- calcNormFactors( dge, method = "TMM")
counts <- cpm(edger, log = TRUE)
selectedFeatures <- rownames( edger )[ apply( counts, 1, function( v ) sum( v >= cpm_value ) ) >= 1/4 * ncol( counts ) ]
# GET HIGH EXPRESSED FEATURES
highExprDge <- dge[ selectedFeatures,, keep.lib.sizes = FALSE ]
# NORMALIZE DATA
normDge <- calcNormFactors( highExprDge, method = "TMM")
tempDge <- normDge
tempDge$counts <- cpm(normDge, log = TRUE)
countDistributionLinePlot(tempDge)
variancePcaPlot(tempDge)
# CREATE DESIGN
design <- model.matrix(eval(parse(text=design_value)), normDge$samples)
design <- reNameDesign(design, normDge)
matrix_value <- strsplit(matrix_value, ",")[[1]]
contr.matrix <- makeContrasts(contrasts = matrix_value, levels=design)
contr.matrix
## Contrasts
## Levels pDC - pDC_unstimulated
## pDC 1
## pDC_unstimulated -1
# PERFORM ANALYSIS
normDge <- estimateDisp(normDge, design)
edfit <- glmFit(normDge, design)
efit <- glmLRT(edfit, contrast = contr.matrix)
efit$DE <- decideTests(efit, p.value = alpha)
summary(efit$DE)
## 1*pDC -1*pDC_unstimulated
## Down 3273
## NotSig 7819
## Up 3130
# CREATE deTab TABLE
deTab <- topTags(efit, n=Inf)$table
deTab <- deTab[order(rownames(deTab)),]
deTab$DE <- c(efit$DE[order(rownames(efit$DE)),])
deTab$adj.P.Val <- p.adjust(deTab$PValue, method = "BH")
deTab <- rename(deTab, "avgLog2CPM" = "logCPM")
deTab <- rename(deTab, "avgLog2FC" = "logFC")
deTab <- rename(deTab, "P.Value" = "PValue")
# CHANGE GENE ID TO SYMBOL IF NESSECARY
# SHOW DGE RESULTS
ma_plot(deTab, NULL)
pValuePlot(deTab)
# SAVE ANALYSIS
normDge$counts <- cpm(normDge, log = TRUE)
save(deTab, normDge, file="analysis.RData")
# INFO
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] reshape2_1.4.3 org.Mm.eg.db_3.10.0 org.Hs.eg.db_3.10.0 AnnotationDbi_1.48.0 igraph_1.2.4.2
## [6] ReactomePA_1.30.0 graphite_1.32.0 DOSE_3.12.0 clusterProfiler_3.14.3 plotly_4.9.1
## [11] ggplot2_3.2.1 broom_0.5.4 scales_1.1.0 tidyr_1.0.2 DESeq2_1.26.0
## [16] edgeR_3.28.0 limma_3.42.1 SummarizedExperiment_1.16.1 DelayedArray_0.12.2 BiocParallel_1.20.1
## [21] matrixStats_0.55.0 Biobase_2.46.0 GenomicRanges_1.38.0 GenomeInfoDb_1.22.0 IRanges_2.20.2
## [26] S4Vectors_0.24.3 BiocGenerics_0.32.0 knitr_1.28 BiocManager_1.30.10 shinyjs_1.1
## [31] shinycssloaders_0.3 shinyWidgets_0.5.0 shinydashboard_0.7.1 shiny_1.4.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.5 Hmisc_4.3-0 fastmatch_1.1-0 plyr_1.8.5 lazyeval_0.2.2 splines_3.6.1 crosstalk_1.0.0
## [8] urltools_1.7.3 digest_0.6.23 htmltools_0.4.0 GOSemSim_2.12.0 viridis_0.5.1 GO.db_3.10.0 magrittr_1.5
## [15] checkmate_1.9.4 memoise_1.1.0 cluster_2.1.0 annotate_1.64.0 graphlayouts_0.5.0 enrichplot_1.6.1 prettyunits_1.1.1
## [22] jpeg_0.1-8.1 colorspace_1.4-1 rappdirs_0.3.1 blob_1.2.1 ggrepel_0.8.1 xfun_0.12 dplyr_0.8.4
## [29] crayon_1.3.4 RCurl_1.98-1.1 jsonlite_1.6.1 graph_1.64.0 genefilter_1.68.0 survival_3.1-8 glue_1.3.1
## [36] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.32.0 XVector_0.26.0 DBI_1.1.0 Rcpp_1.0.3 viridisLite_0.3.0
## [43] xtable_1.8-4 progress_1.2.2 htmlTable_1.13.3 gridGraphics_0.4-1 reactome.db_1.70.0 foreign_0.8-75 bit_1.1-15.1
## [50] europepmc_0.3 Formula_1.2-3 DT_0.11 htmlwidgets_1.5.1 httr_1.4.1 fgsea_1.12.0 RColorBrewer_1.1-2
## [57] acepack_1.4.1 pkgconfig_2.0.3 XML_3.99-0.3 farver_2.0.3 nnet_7.3-12 locfit_1.5-9.1 ggplotify_0.0.4
## [64] tidyselect_1.0.0 rlang_0.4.4 later_1.0.0 munsell_0.5.0 tools_3.6.1 generics_0.0.2 RSQLite_2.2.0
## [71] ggridges_0.5.2 evaluate_0.14 stringr_1.4.0 fastmap_1.0.1 yaml_2.2.1 bit64_0.9-7 tidygraph_1.1.2
## [78] purrr_0.3.3 ggraph_2.0.0 packrat_0.5.0 nlme_3.1-143 mime_0.9 DO.db_2.9 xml2_1.2.2
## [85] compiler_3.6.1 rstudioapi_0.10 png_0.1-7 tibble_2.1.3 tweenr_1.0.1 geneplotter_1.64.0 stringi_1.4.5
## [92] lattice_0.20-38 Matrix_1.2-18 vctrs_0.2.2 pillar_1.4.3 lifecycle_0.1.0 triebeard_0.3.0 data.table_1.12.8
## [99] cowplot_1.0.0 bitops_1.0-6 httpuv_1.5.2 qvalue_2.18.0 R6_2.4.1 latticeExtra_0.6-29 promises_1.1.0
## [106] gridExtra_2.3 MASS_7.3-51.5 assertthat_0.2.1 withr_2.1.2 GenomeInfoDbData_1.2.2 mgcv_1.8-31 hms_0.5.3
## [113] grid_3.6.1 rpart_4.1-15 rmarkdown_2.1 rvcheck_0.1.7 ggforce_0.3.1 base64enc_0.1-3